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Slow crack propagation in ductile, and in certain brittle materials, appears to take place via the 
nucleation of voids ahead of the crack tip due to plastic yields, followed by the coalescence of these 
voids. Post mortem analysis of the resulting fracture surfaces of ductile and brittle materials on the 
/^m-mm and the nm scales respectively, reveals self-affine cracks with anomalous scaling exponent 
C, ~ 0.8 in 3-dimensions and C, ~ 0.65 in 2-dimensions. In this paper we present an analytic theory 
based on the method of iterated conformal maps aimed at modelling the void formation and the 
fracture growth, culminating in estimates of the roughening exponents in 2-dimensions. In the 
simplest realization of the model we allow one void ahead of the crack, and address the robustness 
of the roughening exponent. Next we develop the theory further, to include two voids ahead of the 
crack. This development necessitates generalizing the method of iterated conformal maps to include 
doubly connected regions (maps from the annulus rather than the unit circle). While mathematically 
and numerically feasible, we find that the employment of the stress field as computed from elasticity 
theory becomes questionable when more than one void is explicitly inserted into the material. Thus 
further progress in this line of research calls for improved treatment of the plastic dynamics. 



I. INTRODUCTION 



In this paper we expand on the results of a recent Let- 
ter in which a model was proposed for slow crack prop- 
agation via void formation ahead of the crack due to plas- 
tic yields. Here 'slow' means propagation velocity consid- 
erably smaller than the Rayleigh wave speed. The model 
was motivated by some quantitative studies of fracture 
surfaces, which reveal self-affine rough cracks with two 
scaling regimes: at small length scales (smaller than a 
typical cross-over length ^c) the roughness exponent is 
C « 0.5, whereas at scales larger than the roughness 
exponent is C « 0.8. The second scaling regime is seen to 
have an upper cut-off ^ known as the correlation length. 
Such measurements were reported first for ductile mate- 
rials (like metals) where is of the order of 1 fim. 0, Q , 
and more recently for brittle materials like glass, but with 
a much smaller value of ^c, of about 1 nm Similar ex- 
periments conducted on 2-dimensional samples reported 
rough cracks with large-scale exponents C, ~ 0.65 ± 0.04 
lilEH- The exponent C « 0.5 is characteristic of uncor- 
related random surfaces, but higher exponents indicate 
the existence of positive correlations naturally, the 
experimental discovery of such correlated "anomalous" 
exponents attracted considerable interest with repeated 
attempts to derive them theoretically 9, 10, nj . Ref . jj 
presented a quantitative model for self-affine fracture sur- 
faces based on elasticity theory supplemented with con- 
siderations of plastic deformations. Focussing on infinite 
2-dimensional materials, Ref. Tl followed the qualitative 
picture presented recently in |12| , see Fig. ^ In this pic- 
ture there exists a "process zone" in front of the crack tip 
in which plastic yield is accompanied by the evolution of 
damage cavities. A crucial aspect of this picture is the ex- 
istence of a typical scale, ^c, which is roughly the distance 
between the crack tip and the first void, at the time of 
the nucleation of the latter. The voids are nucleated un- 
der the influence of the stress field (Jij{r) adjacent to the 
tip, but not at the tip, due to the existence of the plastic 



(^ = 0.5 roughness of the cavities , ^ 
(,=0.65 roughness of the structure formed by damage cavities 



FIG. 1: The fracture scenario suggested in [l^ . This scenario 
had been documented in detail in corrosive glass fracture, and 
also more recently in the fracture of paper plj |. 



zone that cuts ofi^ the purely linear-elastic (unphysical) 
crack-tip singularities. The crack grows by coalescing the 
voids with the tip, creating a new stress field which in- 
duces the nucleation of new voids. In the picture of 
the scale is also identified with the typical size of the 
voids at coalescence. A consequence of this picture is 
that the roughening exponent C w 0.5 corresponds to the 
surface structure of individual voids, whereas the large- 
scale anomalous exponent has to do with the correlation 
between the positions of different voids that coalesce to 
constitute the evolving crack. Ref. provided first a 
theory for the scale and second, demonstrated that 
the positions of consecutive voids are positively corre- 
lated. In this paper we amplify on these results. 

In Sect, m we discuss again the model of Ref. P|, 
expand the theoretical presentation, and test the robust- 
ness of the roughening exponent that this model predicts. 
Since the roughening exponent found is larger than 0.5, 
this indicates the existence of positive correlations be- 
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tween crack increments. To illuminate these long range 
correlation we demonstrate in Sect. 11111 that the analytic 
structure of the theory dictates the existence of power- 
law correlations between height fluctuations in the crack 
and the value of the mode II stress intensity factor K^-^ 
at the tip. In Sect. |EI we construct a model with two 
voids ahead of the crack. After setting up the problem, 
we describe in some detail the mathematical apparatus 
that employs conformal maps from the annulus to the 
doubly connected region of a crack with a void which 
allows the computation of the stress field around such 
a configuration. In Sect. we present results of the 
two-void model, and discuss the relevance of the results 
to the growth and roughening of cracks. The conclu- 
sion is that since the details of plastic deformations are 
not well understood the physics of crack growth is bet- 
ter described by the one-void model than the two-void 
model; the stress field computed for a crack with one 
void ahead is physically acceptable as long as elasticity 
theory is relevant, but when the first void appears due to 
plastic events a correct determination of the stress field 
should include a better handling of the plastic zone. This 
must await future improvement of our understanding of 
plastic dynamics. 



II. ONE- VOID MODEL 

A. The plastic zone and void nucleation 

A simple model for can be developed by assuming 
the process zone to be properly described by the Huber- 
von Mises plasticity theory This theory focuses on 



the deviatoric stress = try - 



^TicrSij and on its invari- 
1 , 



ants. The second invariant, J2 = ^SijSij, corresponds to 
the distortional energy. The material yields as the dis- 
tortional energy exceeds a material-dependent threshold 
a^. The fact that we treat this threshold as a constant, 
independent of the state of deformation and its history, 
implies that we specialize for "perfect" plasticity. In 2 
dimensions this yield condition reads 
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Here ai^2 are the principal stresses given by 
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In the purely linear-elastic solution the crack-tip region 
is where high stresses are concentrated (in fact diverging 
near a sharp tip). Perfect plasticity implies on the one 
hand that the tip is blunted, and on the other hand that 
inside the plastic zone the Huber-von Mises criterion 
is satisfied. The outer boundary of the plastic zone will 
be called below the "yield curve", and in polar coordi- 
nates around the crack tip will be denoted R{6). 



Whatever is the actual shape of the blunted tip its 
boundary cannot support normal components of the 
stress. Together with Eq. J^l this implies that on the 
crack interface 
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On the other hand, the linear-elastic solution, which is 
still valid outside the plastic zone, imposes the outer 
boundary conditions on the yield curve. Below we will 
compute the outer stress field exactly for an arbitrarily 
shaped crack using the recently developed method of iter- 
ated conformal mappings |l5j | . For the present argument 
we will take the outer stress field to conform with the 
universal linear-elastic stress field for mode I symmetry. 



Gij (r, B) = 



(4) 



For a crack of length L with cr°° being the tensile load at 
infinity, the stress intensity factor is expected to scale 
like ^ a°"^/Tj. Using this field we can find the yield 
curve R{9). Typical yield curves for straight and curved 
cracks are shown in the insets of Figs. |31andEl 

The typical scale follows from the physics of the nu- 
cleation process. It is physically plausible that void for- 
mation is more susceptible to the growth of hydrostatic 
tension than to distortional stresses. We assume that 
void nucleation occurs where the hydrostatic tension P, 
P = ^Trcr, exceeds some threshold value Pc- The hydro- 
static tension increases when we go away from the tip and 
reaches a maximum near the yield curve. To see this note 
that on the crack surface P = (cf. Eq. Q). On 

the yield curve we use Eq. Q and the Huber-von Mises 
criterion together to solve the angular dependence of the 
hydrostatic tension in units of . It attains a maximal 

value of VScTy and is considerably higher than for 
a wide range of angles. On the other hand the linear- 
elastic solution Q implies a monotonically decreasing P 
outside the yield curve. We thus expect P to attain its 
maximum value near the yield curve. This conclusion is 
fully supported by finite elements method calculations, 
cf. 16]. Finally, since the nucleation occurs when P ex- 
ceeds a threshold Pc, this threshold is between the limit 
values found above, i.e. ^a.^<Pc<VScr^ . The void will 
thus appear at a typical distance ^c, see Fig. |21 An im- 
mediate consequence of the above discussion is that is 
related to the crack length via: 



a2 



(5) 



Note that is not a newly found length scale; it is 
the well known scale of the plastic zone 17]. Its iden- 
tification with the cross-over length between two scaling 
behaviors of the crack roughening is however new. This 
stems from the proposition that positive correlations ap- 
pear only between the positions of nucleated voids. Be- 
low one enters the regime of plastic processes whose 
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FIG. 2: A forward direction profile of tfie fiydrostatic tension 
P in units of a^. On the craclt P — V3/2 and it attains 
a maximum of ^/3 on the yield curve. The threshold line 
indicates a value of Pc such that ^ < Pc< \/3. The typi- 
cal length is shown. Other directions exhibit qualitatively 
similar profiles. 

theory is far from being settled. We should also comment 
that it is possible that positive correlations appear even 
below the scale of the plastic zone since experiments in- 
dicate that several voids nucleate within the plastic zone 
[T^IT3| . We should therefore consider the estimate in Eq. 
Q as an upper bound on Finally, the fact that the 
plastic zone size scales with as proposed in Eq. Q 
results from the assumption of perfect plasticity, i.e. that 
cTy is independent of the state of deformation and its his- 
tory. This is not true for real materials; usually is not 
sharply defined; it can increase with plastic deformations 
|l4j . This phenomenon, known as "work-hardening" or 
"strain-hardening" , might introduce other dependencies 
on and include other length scales that are related to 
the plastic deformations. We do not take such issues into 
account in this simple model. 

Naturally, the precise location of the nucleating void 
will experience a high degree of stochasticity due to 
material inhomogeneities. Since we do not know from 
first principles the probability distribution for void 
formation, we consider in our model below two possible 
distribution functions. In all cases nucleation cannot oc- 
cur if P < Pc- For P>Pc the void occurs with probability 



sponsible later for the roughening of the crack. For com- 
parison examine also the pdf 's for a general crack which 
are shown in Fig. |H| There the symmetry is lost: correla- 
tion to previous steps create a preference for the upward 
direction. This source of positive correlations is discussed 
below in greater detail. 
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FIG. 3: Panel (i): the tip of a straight crack and the yield 
curve in front of it. 

Panel (ii): three probability distribution functions calculated 
for the configuration in (i). The abscissa is 0, the angle mea- 
sured from the crack tip as seen in panel (i). The ordinate 
is the normalized probability (per unit 6) to grow in the 9 
direction. The distributions are symmetric and wide enough 
to allow deviations from the forward direction. For all the 
curves = 6. For curve (a) p{9) oc exp[(P — Pc)] — 1 and 
^ = 8, for curve (b) p{9) oc exp[0.2(P-Pc)] - 1 and ^ = 6 
and for curve (c) p(^) (x P — Pc and -p^ = 6. 



B. Crack Propagation 



P GC P-Pc , (6) 

P cx exp[a(P - Pc)] - 1 . (7) 

In the exponential case we considered two different val- 
ues of a. In Fig. O we show three such pdf 's as they 
appear for a perfectly straight crack. We note that these 
distributions are symmetric about the forward direction. 
Nevertheless they have sufficient width to allow devia- 
tions from forward growth. These deviations will be re- 



Each growth step in our model is composed of two 
events. Firstly the material yields near the crack tip, 
creating a plastic zone with a void growing somewhere 
at the zone boundary. Secondly the crack tip and the 
void coalesce. We note that there is a separation of time 
scales between these two events. The first is slow enough 
to be governed by a quasi-static stress field. The second 
event occurs on a shorter time scale. It is clear that we 
forsake in the one-void model any detailed description 
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of the geometry on scales smaller than ^c- Any relevant 
scaling exponent that will be found in this model will 
refer to roughening on length scales larger than ^c- In 
experiments it appears that several voids may nucleate 
before the coalescence occurs, and in the next section we 
will explore models with two voids ahead of the crack. 
In the one-void model, the physical process in which the 
crack coalesces with the multiple voids ahead of it is sub- 
stituted by a single void coalescence with the crack. 

In spite of the simplification in dealing with only one 
void per step, it was demonstrated in J,] that the one- 
void model induces positive correlations between consec- 
utive void nucleations, leading eventually to an anoma- 
lous roughness exponent larger than 0.5. Clearly, even 
this simple model requires strong tools to compute the 
stress field around an arbitrarily shaped crack, to deter- 
mine at each stage of growth the location of the yield 
curve and nucleating randomly the next void according 
to the probability distributions discussed above. In a re- 
cent work we have developed precisely the necessary tool 
in the form of the method of iterated conformal mappings 

m 

In the method of iterated conformal mappings one 
starts with a crack for which the conformal map from 
the exterior of the unit circle to the exterior of the crack 
is known. (Below we start with a long crack, in the form 
of a mathematical branch-cut of length 10000, and is 
of 0(10)). We can then grow the crack by little steps in 
desired directions, computing at all times the conformal 
map from the exterior of the unit circle to the exterior of 
the resulting crack. Having the conformal map makes the 
exact calculation of the stress field (for arbitrary loads at 
infinity) straightforward in principle and highly afford- 
able in practice. The details of the method and its ma- 
chine implementations are described in full detail in 'l5] . 
In the next section we present the theory in great detail 
for the two-void model, and avoid the repetition here. 
We should just stress that the method naturally grows 
cracks with tips of finite curvature, and each step adds 
on a small addition to the tip, also of a finite size that is 
controlled in the algorithm. 

Having the stress field around the crack we can readily 
find the yield curve and the physical region in its vicinity 
where a void can be nucleated. Choosing with any one 
of the probability distributions described above, we use 
this site as a pointer that directs the crack tip. We then 
use the method of iterated conformal mappings to make 
a growth step to coalesce the tip with the void. Natu- 
rally the step sizes are of the order of ^c- In Fig. 01 we 
present the actual step sizes as computed with the pdf 
©, as a function of the crack length. The linearity (in 
the mean) in L is obvious. Note that the fluctuations 
about the mean are strongly dependent on the pdf, and 
could in principle be used to experimentally deduce the 
'correct' pdf by reverse engineering. We reiterate that 
this model forsakes the details of the void structure and 
all the length scales below ^c- Since we are making linear 
steps below we anticipate having an artificial scal- 



ing exponent C = 1 for scales smaller than fc- This is 
clearly acceptable as long as we are mainly interested in 
the scaling properties on scales larger than ^c- 
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crack length (units of bump radius) 

FIG. 4: The step size obtained using pdf l|SJ (see section 
rm . vs. the crack length. Note that these results agree with 
Eq.@ i.e. the step size grows linearly with the crack length. 
The units here are the "bump" radius which is introduced 
explicitly in Eq. lUTl . 

In Fig. [Slwe present two typical cracks that were grown 
using this method. Both cracks were initiated from a 
straight crack of length 10000. The upper crack was 
grown using the broader exponential pdf of Fig. |5| curve 
(b). The lower crack was grown with the narrower pdf 
of Fig. El curve (a). Clearly, the upper crack exhibits 
stronger height fluctuations, as can be expected from the 
wider pdf and the choice of parameters. For the lower 
crack forward growth is much more preferred. In the 
upper crack the positive correlations between successive 
void nucleation and coalescence events can be seen even 
with the naked eye. This is precisely the property that 
we were after. A neat way to see this tendency is in 
the pdf 's as they are computed on the yields crack for a 
typical, rather than straight, crack. In Fig. [Slwe show 
these pdf's for the crack whose yield curve is shown in 
the upper panel. We see that now the symmetry of the 
pdf's is lost, and positive values of 9 are preferred. This 
is the source of positive correlations that eventually give 
rise to the anomalous roughening exponent. This is born 
out by the measurements of the scaling exponent that we 
discuss next. 

A quantitative measurement of the positive correla- 
tions is the roughening exponent, that we compute as 
follows. Measuring the height fluctuations y{x) in the 
graph of the crack, one defines h{r) according to 

hir) EE (Maxivix)}^^^^.^^^ - Mzn {y(i)}^<,<^^,,)^ . 

(8) 

For self-affine graphs the scaling exponent ^ is defined 
via the scaling relation 

h{r) - . (9) 
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FIG. 5: Two typical cracks generated with our model. Note 
the different scales of the abscissa and ordinate, and that the 
lower crack had been translated by -300. The upper crack 
exhibits two decades of self-afiine scaling with a Hurst expo- 
nent 0.64. The lower crack has smaller standard deviation 
and therefore a shorter scaling range. Nevertheless it appears 
that in its shorter scaling range it exhibits an exponent that 
is very close to the upper crack. 



In Fig. [71 we present a typical log-log plot of h{r) vs. 
r, in this case for the two cracks in Fig. [S] with power- 
law fits of C = 0.64 and 0.68 respectively. Indeed as 
anticipated from the visual observation of Fig. |Slthe ex- 
ponent is higher than 0.5. It turned out that all the 
cracks grown by our algorithm gave rise to scaling plots 
in which a scaling range with ( = 0.66 ± 0.03 is clearly 
seen. When the pdf allowed for a sizeable standard de- 
viation, the cracks gave a very nice scaling plot with at 
least two decades of clear anomalous scaling. When the 
standard deviation was small, the scaling range was more 
meager, as seen in Fig. |7| It is interesting to stress that 
the anomalous scaling exponent appears insensitive to 
the pdf used (although the extent of the scaling range 
clearly depended on the pdf). We note that our mea- 
sured scaling exponents are very close to the exponents 
observed in 2-diniensional experiments. (Of course we 
cannot expect a 2-dimensional theory to agree with 3- 
dimensional experiments - the scaling exponents are, as 
always, dimension-dependent). In addition the value of 

does not effect the scaling properties of a crack, i.e. it 
doesn't seem to matter how long the step is, so long as a 
wide distribution of angles is allowed. 

Growing directly at the tip of the crack results in a very 
strong preference for the forward direction, meaning that 
a step up will most likely be followed by a step down, 
and vice versa, as shown in [T^ . The introduction of 
the physics of the plastic zone results in creating a finite 
distance away from the tip to realize the next growth 
step. Another crucial issue is the existence of long range 
correlations. Since this aspect was not made clear so far, 
we turn now to a discussion of the origin of power law 
correlations. 
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FIG. 6: Panel (i): the tip of a "rough" crack and the yield 
curve in front of it. Panel (ii): three probability distribution 
functions calculated for the configuration in panel (i). The 
abscissa is 9, the angle measured from the crack tip as seen in 
panel (i) . The ordinate is the normalized probability (per unit 
0) to grow in the 6 direction. The pdf 's are those used in Fig. 
121 using the same parameters. Note the upward preference in 
all the pdf 's due to the broken symmetry. 



III. LONG RANGE POSITIVE CORRELATIONS 
IN FRACTURE 



The fact that the cracks generated by our model ap- 
pear self-affine with Hurst exponent C > 0.5 implies that 
the physical mechanism underling the crack growth is a 
long range positive correlation process. We can gain in- 
tuition about the origin of the long range correlations by 
employing some known analytic results. Consider a long 
mode I straight crack spanning the interval [—L, 0]. Sup- 
pose now that the crack shape is perturbed by a small 
out of plane fluctuation of the form eip^x), where e > 
is small. In the presence of the perturbation the crack 
attains a small shear component K^^ 7^ at its tip. A 
first order perturbation analysis in the amplitude e re- 
veals that 19] 
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FIG. 7: Calculation of the anomalous roughening exponent. 
The slopes of the dotted lines are 0.64 for the upper plot 
(curve a) and 0.68 for the lower (curve b). Note that the initial 
scaling with slope 1 is relevant for length scales smaller than 
^c- This scaling is unphysical, resulting from our algorithm 
that connects the crack tip to a void by a straight line. 



where the superscript in a^J{x,0) refers to the solution 
in the absence of a perturbation. Note that due to the 
fact that the crack is long we could set the lower limit 
of integration to — oo. Let us consider a positive pertur- 
bation i^i^) that is symmetric around —r and decays to 
zero at a typical distance 6 from — r. Since in our config- 
uration cri?(a;, 0) — —a°° for a; < 0, it is straight forward 
to show that Eq. ifTUIl yields 



r3/2 



(11) 



for r ^ (5. Note that K^^ is negative. In order to under- 
stand the effect of the perturbation on the probability 
distribution function for the next void nucleation in our 
model, we recall that this probability is determined by 
the hydrostatic tension P = ^Trcr. In the case of a pure 
tensile stress, K^^ = 0, P is symmetric around 9 = 
and the probabilities of nucleating a void at positive or 
negative angles are the same. In the presence of a small 
negative shear component, K^^ < 0, this picture changes. 
The maximal hydrostatic tension is obtained at 



(12) 



Since to first order in e the mode I stress intensity factor 
if J ~ a°°y/lj is unchanged obtain that the peak 

of the distribution is shifted from zero to 
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^1/2^3/2 
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(13) 



This relation shows that as a result of a positive pertur- 
bation in the crack shape at a distance r behind its tip the 



probability to nucleate a void at a positive angle relative 
to the forward direction is higher than the probability 
to nucleate a void at a negative angle. Moreover, this 
positive correlation is long ranged, decaying as r~^^^. 

The presence of long range correlations is reassuring, 
since they are a must for the existence of a roughening 
exponent larger than 0.5. Note however that the above 
result does not determine in any direct way the numerical 
value of the roughening exponent itself. The actual ex- 
ponent results from the cumulative effect of many height 
fluctuations, and at present we do not have an analytic 
theory predicting the numerical value of this exponent. 
Contrary to self similar fractal growth patterns, where 
the fractal dimension D can be computed from the knowl- 
edge of the first Laurent coefficient of the conformal map 
|20|, in self afhne graphs it is not obvious how to extract 
the roughening exponent from the properties of the con- 
formal map. At present we are bound to the laborious 
process of actually growing the crack and measuring the 
exponent. Needless to say this is theoretically unsatisfac- 
tory, and new ideas on this issue should be very welcome. 



IV. TWO- VOID MODELS: MATHEMATICAL 
FORMULATION 

In this section we address the physical process of nu- 
cleation of a second void in front of the crack tip. To 
improve upon the one-void model we need to develop 
techniques to compute the exact stress field around a 
crack with one void ahead (or a crack with two voids 
ahead, etc.). Clearly, methods based on conformal maps 
from singly connected regions (i.e. the unit circle) can- 
not suffice for this purpose. Since conformal map tech- 
niques from doubly or multiply connected regions are far 
less familiar, and since the computation of stress field 
around doubly connected regions is interesting by itself, 
we present the necessary techniques in some detail. Con- 
siderable attention will be paid to the accuracy and the 
efficiency of the calculation. In subsection II V Bl we re- 
view the basis of the relevance of conformal maps to the 
solution of the bi-Laplace equation. This goes back to 
Mushkelishvili's series-expansion method (2l|. We ex- 
tend Mushkelishvili's method to doubly connected re- 
gions, using conformal maps from the annulus to the re- 
quired doubly connected region. In Appendix^we elab- 
orate the case of two circular holes. This case is solvable 
analytically using bipolar coordinates , and therefore 
provides a unique testing ground for the precision of our 
method in a limiting case. 

In this section we solve for the stress field in an un- 
bounded planar doubly connected region. Such calcula- 
tions exist in the literature; one general method is known 
as the Schwarz alternating method. In this method one 
solves the simply connected problems successively and 
superimposes them in such a way that the boundary con- 
ditions are satisfied when the procedure has converged 
(see for example 2§l). We will assess our method by com- 
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paring with an analytical solution for two equally sized 
circular holes using bipolar coordinates ^ . Our method 
is quite genearal, allowing us to solve for the stress field 
in any doubly connected region by explicitly solving the 
elastic equations in a doubly connected geometry. The 
advantages of our method arc that it allows freedom in 
choosing the shapes of the boundaries and it allows cal- 
culation of the stress field near highly singular shapes 
such as long cracks. The method can be implemented to 
very high precision as will be shown below. 



A. Equilibrium equations for the stress field in a 
doubly connected infinite medium. 

The theory of elastostatic fracture mechanics in brittle 
continuous media is based on the equilibrium equations 
for an isotropic elastic body psf 



0. 



(14) 



For in-plane modes of fractures, i.e. under plane-stress 
or plane- strain conditions, one introduces the Airy stress 
potential U{x,y) such that 



dy 



2 '"^2/ 



dxdy 



; CTyy 



(15) 



Thus the set of Eq. H14|) . after simple manipulations, 
translate to a Bi-Laplace equation for the Airy stress po- 
tential U{x,y) ;25j 



(16) 



with the prescribed boundary conditions on the crack 
and on the external boundaries of the material. At this 
point we choose to focus on the case of uniform remote 
loadings and traction-free crack boundaries. This choice, 
although not the most general, is of great interest and 
will serve to elucidate our method. Other solutions may 
be obtained by superposition. Thus, the boundary con- 
ditions at infinity, for the two in-plane symmetry modes 
of fracture, are presented as 

(Jxxioo) = ■,(Tyy{oo) = a°° ; axy (oo) = Mode I (17) 
cTxxioo) = ■,ayy{oo) = ;(Ta;y(oo) = a°° Mode II . 

In addition, the free boundary conditions on both bound- 
aries (of crack and void) are expressed as 



O'xn(s) = CTynis) = 



(18) 



where s is the arc-length parametrization of the bound- 
aries and the subscript n denotes the out-ward normal 
direction. 

The solution of the Bi-Laplace equation can be written 
in terms of two analytic functions 0(z) and r]{z) as 



C/(x,2/)=5R[z^(z) + r;(z)] 



(19) 



In terms of these two analytic functions, using Eq. I|15|l . 
the stress components are given by 

ayy{x,y) = n[2^'{z) + z^"{z) + 7j"{z)] 
cTxx{x,y) = ^[2^\z)~zip"{z)-r^"{z)] 
a^y{x,y) = (20) 

In order to compute the full stress field one should first 
formulate the boundary conditions in terms of the ana- 
lytic functions ip{z) and ri{z). The boundary conditions 
Eq. ((TH)l can be rewritten, using Eq. ifTCjl . as plf 



dt 



dx dy 







(21) 



Where dt is the tangential derivative along the bound- 
aries. Condition (|21|l must hold for each of the two 
boundaries separately. Note that we do not have enough 
boundary conditions to determine U{x, y) uniquely. In 
fact we can allow in Eq. (|19|l arbitrary transformations 
of the form 



ip + iCz + 7 



tp + J , ^ = 7]' 



(22) 



where C is a real constant and 7 and 7 are complex 
constants. This provides five degrees of freedom in the 
definition of the Airy potential. It is important to stress 
that whatever the choice of the five freedoms, the stress 
tensor is unaffected; see for an exhaustive discussion 
of this point. We will explain below how to take advan- 
tage of these freedoms to make the formulation simpler. 
When the domain is doubly connected, the traction-free 
conditions H21() can be written using H19|l as 



(p{z) + zip'{z) + 'ip{z) = Di for z eCi 



(23) 



and 



^{z) + z^'{z) + ij{z)^D2 for zeC2 (24) 

where Ci and C2 are the two boundary curves, and Di, 
D2 are complex constants that are eventually uniquely 
determined by the solution of the problem. To proceed 
we represent ip{z) and i^{z) in Laurent expansion form. 
Note that (p{z) and ip{z) have a poles inside both the 
boundaries and therefore do not have a Laurent expan- 
sion around infinity which is valid everywhere in the com- 
plex plane. But for \z\ > R where i? is a radius that 
encloses both the boundaries the following expansion is 
valid since f{z) and iJj{z) have no poles in this region. 



Ifiiz) = ipiZ + ipQ + ip_i/z + ip_2/z^ + 



(25) 



This form is in agreement with the boundary conditions 
at infinity that disallow higher order terms in z. Using 
the boundary conditions ((T7|l . we find 



^1 



; 



01 



Mode I , 
Mode II . 



(26) 
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Where one of the freedoms in H22(l was used to choose tpi 
to be real, using the real constant C in H22I) . The four 
remaining freedoms will allow us later on to fix lpq and 
^0 in a convenient way. 



B. Application of conformal maps. 

In order to enable the calculation of the stress field 
around an arbitrarily shaped crack and a void, we con- 
formally map the annulus (having its outer radius set to 
one, i.e. the annulus is p < r < 1) onto the required 
doubly connected domain, see Fig. |S1 A well known 
fact is that simply connected domains can be conformally 
mapped to any other simply connected domains (relying 
on the Riemann mapping theorem). However, when deal- 
ing with doubly connected domains there is an invariant 
quantity, called the modulus (sometimes recasted as the 
extremal distance), which is preserved under conformal 
mappings. As a result only doubly-connected domains 
with the same modulus can be connected via a confor- 
mal map. For an annulus the conformal modulus is just 
the ratio of the inner radius and the outer radius, so that 
for the p < r < 1 annulus, the modulus is simply p. For 
that reason the specific annulus which is taken as the 
domain to be mapped onto the required crack -f void do- 
main cannot be just any annulus, but has to be chosen 
correctly. 
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FIG. 8: Illustration of how the conformal map $(0;) operates. 
The unit circle is mapped to Ci and the inner circle to C2. 
The point a is mapped to infinity. 

Suppose that we have such a conformal map $(tj) (ex- 
amples are presented in the following subsections) that 
maps the w-annulus domain onto the required physi- 
cal z-plane with a crack and a void in front of it. We 
employ this map to find the solution of the stress field 
in the physical domain (we also introduce the notation 
x[z) = $~^(z) for the inverse map). Due to Eq. (fT^ . 
knowing the solution to the Bi-Laplace equation in the 
annulus does not immediately provide the solution as the 
Bi-Laplace equation in the physical domain through a 
simple application of the conformal map, since in contrast 
to the Laplace equation, the Bi-Laplace equation is not 
conformally invariant. Nevertheless, the conformal map- 
ping method can be extended to non-Laplacian problems 
and provides a clear simplification of the problem since 
the boundary conditions are much easier to impose on 
a circular boundary than on the physical boundary. We 
begin by writing our unknown functions ip{z) and ip{z) 



in terms of the conformal map 

'P (z) = If ix (z)) and ijj (z) = -ip (x (z)) (27) 

Now, (fiiuj) and ip{uj) are just analytical functions in the 
annulus, apart from a simple pole located at a point 
which is mapped to infinity in the z-plane, cf. Fig. |S1 
Actually, a closer inspection leads us to use the fact that 
in the z-plane ip{z) and i'iz) have truncated expansions 
as given in Eqs. H25|l . Thus, we expect ip{uj) and Tpiuj) 
to be of the general form 



n— — 00 



and 



n— — 00 



where for Mode I fracture we can identify 

A = Lpi = — ^ and B — th-i = — ^ 



(28) 



(29) 



(30) 



just like in Eq. (|26|l (from here on we use only Mode I 
boundary conditions at infinity). In contrast to simply 
connected domains, here we have positive as well as neg- 
ative powers of lu in the expansion. At this point we use 
the four remaining freedoms (7 and 7 in H22I) ') to choose 
ipo and ipo such that 



ipo = and ipQ — 



(31) 



To impose boundary conditions on the outer (unit) circle 
of the annulus (i.e. |w| — 1) we write Eq. H23|) in the uj 
plane. This yields (using e — e^^) 



oc 



(32) 



11— —00 



n— — 00 



Imposing boundary conditions on the inner circle of the 
annulus (i.e. \uj\ = p) we write Eq. H24|l in the uj plane. 
This gives (using uj = pe = pe'^) 



[$(p£) + $(pe) 



E 



fnP £ + 



(33) 



n— — 00 



^^P'^ E nrnp-'e-+'+ E 



= D2. 



n— — 00 



{pe 

To proceed we must Fourier Transform the functions 
^{ijj) and <I>(tj)/$'(ijj) on the boundaries of the annulus. 
We use the following notation: 



$(e) 



E 



out n 



9 



4>(e) ^^j^j ^ indices with 'in' correspond to the inner circle of radius p. 

~ 2-^ " ' Note that the Fourier coefficients on the inner boundary 

are not the same as those on the outer boundary, this is 
because 3>(w) has a pole inside the annulus. Also note 
that contrary to the singly connected case, the expansion 
of ^{uj) goes to infinity in both positive and negative di- 
^(pe) _ ^m/ rr^^\ rcctious. Inserting Eq. into Eos. Ipl^ and (p^ . and 

$' (Qf] ^ n \y ) ■ \ J gathering powers of e we obtain the following infinite set 

of equations for the coefficients (fn , ipn and the unknown 



$(P6) = V C(pe)" 



- — oc 



Here indices with 'out' refer to the outer unit circle, and constants Di and D 



2- 



(pn + + ^K^+n-lWk = (C* + C°4J) + -^n.O^l , n = . . . , OO , (35) 

A;— — OO 

^ ^oo . , 

ipnP +1p-nP + 2^ kb^.^„_-i^p ^ ipk = 2~ V " ^ +C-n/0 j + d„,o£'2 , ^ = -OO , . . . , OO . 



fc— — OO 



This set of equations is well-posed and can be solved 
in various ways. The simplest method is a truncation 
scheme in which one neglects higher order terms in Eqs. 
(|28|l and (|29|l . (taking only a finite subset of coefficients 
ifn and -(/;„, n 7^ 0) and just enough equations for solving 
those coefficients (as well as Di and I?2)- When more 
and more coefficients (/?„ and ipn are taken this scheme 
converges to the exact solution. The efficiency and rate 
of convergence of this simple scheme will be examined 
below. 

The calculation of the Laurent expansion form of ip{uj) 
and ip{i-o) provides the solution of the problem in the uj- 
plane. Still, one should express the derivatives of (p{z) 
and ri{z) in terms of <^(w) and ipi'-^) and the inverse map 
xiz) to obtain the solution in the physical z-planc. This 
is straightforward: 

^'(z) = 0'[x{z)] x'{z) 

^"{z) = ^"[x{z)] [x'{z)]' + ^'[x{z)]x"{z) 
il"{z) = ^'{z) = ^'[x{z)] x'{z). (36) 

Upon substituting these relations into Eq. (|20|l one can 
calculate the full stress field for an arbitrary doubly con- 
nected infinite region. 

In Appendix we present the application of this for- 
malism to the case of two circular holes where the con- 
formal map $([j) and its fourier coefficients are known 
analytically. In this case one achieves extremely accu- 
rate solutions that can be used as testing grounds for the 
truncation method that is always available even when the 
fourier coefficients of the conformal map are not given 
analytically. 



C. Conformal map for arbitrary crack with circular 
void ahead 

Our next task is to find the conformal map ^(lu) from 
the interior of the annulus to the exterior of an arbitrarily 
shaped crack and a void near its tip. We construct this 
conformal map by composing three auxiliary maps. The 
properties of the stress maps are as follows: 

1. Auxiliary Map 1 

Consider the map 4>i{^) given by, 

, , auj — 1 , ^ 

M^) = , 37) 

Lo — a 

a e 3? and < a < 1 . 

(/ii(cj) maps the annulus, i.e. p < |cj| < 1, onto the exte- 
rior of the unit circle and an additional circle on the right 
as is exemplified in Fig. |5| Note that the unit circle is 
mapped onto itself and the inner circle is mapped onto 
the circle on the right. The map (pi has two parameters, 
a and p, the first appearing in its definition (the point 
mapped to infinity), and the other in its domain of def- 
inition. Both a and p are determined by the radius and 
the center of the rightmost circle (see Fig. O) as follows: 

^ _ Xl + X2 

1 + XlX2 + ^{x{ - \){xi - 1)' 
-\^XiX2 - ^J{x\ - \){xl - 1) 

p — , (dSj 

X\ — X2 

where xi and X2 are defined in Fig. O Note that the 
following inequality must hold: 

< p < a < 1 . (39) 
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(0 plane 



FIG. 9: Illustration of how the conformal map (j>i operates. 

maps the annulus into the exterior of the unit circle and 
a circle on its right. The unit circle boundary is mapped to 
itself and the inner circle of radius p is mapped to the circle 
on the right with radius , 






The inverse mapping is given by 

= ^ . (40) 

z ~ a 

Note that the inverse mapping is exactly the same symbol 
as the direct one. 



FIG. 10: Illustration of how the conformal map ^{lu) oper- 
ates. 01 maps the annulus onto the exterior of two circles, (1)2 
rotates the whole plane and (^3 maps the left circle into the 
desired crack shape, leaving the right circle almost unchanged 
in its shape. 



2. Auxiliary Map 2 

The second map 4>2 (i^) is given by, 

(l)2{uj) = Luexp{i6) . (41) 

(j>2{^^) rotates the plane by an angle 6 relative to the real 
axis. The inverse mapping is given by 

(/)~i(z) = zexp{-ie) . (42) 

3. Auxiliary Map 3 

The role of ^3 is to map the exterior of the unit circle to 
the exterior of an arbitrary crack shape. Assume for now 
that we have such a map at hand; in subsection IIV C 51 
we present the explicit derivation of this map using the 
tools of iterated conformal map. At this point consider 
the composition of all three maps. 



annulus into the exterior of two circles; then (j)2 is applied 
to allow the circle on the right to be rotated with respect 
to origin. Finally, (f)^ is applied to map the exterior of 
the unit circle to the exterior of an arbitrary crack shape. 
In total, ^{lu) maps the interior of the annulus to the 
exterior of a crack and void, such that the outer boundary 
of the annulus is mapped to the boundary of the crack 
and the inner boundary is mapped to the void boundary. 
Notice that since (f>3 acts on the whole plane, it affects the 
void ,i.e. the circle on the right in Fig. ^| as well as the 
crack shape. An important property of the mapping we 
suggest is that for all the configurations we are interested 
in, applying ^3 does not change the shape of the void in 
an appreciable way, i.e. the void remains almost circular. 

In order to create the mapping for a given crack and 
void configuration one needs a set of points describing 
the crack's path and the void's radius R and center zq. 
First one constructs 03 according to the desired crack 
shape (see section llV C 5|) . What is left is to obtain the 
values of a, p and 9 (see auxiliary maps 1 & 2). First we 
find xi, X2 and 9, 



4. Composition of the basic maps 



arg(,</)3 



(03"' (^0)) 



(45) 



The desired mapping ^(w) from the annulus to the 
exterior of a crack and void is given by. 



$(w) = 03('/>2(0l(w))) ■ 

The inverse map x{z) is given by. 



(43) 



X{z) = ^-\z) ^ ^^\<j^^\^^\z))) . (44) 

The composition of the three auxiliary maps is illustrated 
in Fig. 1101 First 01 is applied and maps the interior of the 



Xi + X2 

2 

X\ - X2 



= I(03"'(^O))| , 

= |(0,^i(zo + i?)-03"'M)|. (46) 



One can verify that using the above values for xi, X2 and 
0, a void with radius R and center zq is obtained in the 
z plane. Substituting x\ and X2 in Eq. (|38|l we obtain 
the values of p and a. 
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FIG. 11: Example of how to construct the conformal mapping 
along a line. 



5. Conformal Map from the unit circle to an arbitrary 
crack shape 

In this subsection we complete the identification of the 
auxiliary map 03 . For the purpose of being reasonably 
self-contained we reiterate here some aspects of the ma- 
chinery of conformal maps. The essential building block 
in the present application, as in all the applications of 
the method of iterated conformal maps is the fundamen- 
tal map 4>\^0 that maps the exterior circle onto the unit 
circle with a semi-circular bump of linear size \/X which 
is centered at the point e'^. This map reads (2^ : 



2w 



1 + w + w [ 1 + 



1 2 1- A 



w'^ w 1 + X 



1/2 



The inverse mapping (PqIq is of the form 



\z - VTTA(z2 - 1) 
1 - (1 + A)z2 - 



(47) 

1/2 

(48) 
(49) 



By composing this map with itself n times with a judi- 
cious choice of series {Ok\k=i and {\k}n=i '^ill con- 
struct 

4>(")('^) that wiU map the exterior of the circle 
to the exterior of an arbitrary simply connected shape. 
To understand how to choose the two series {Ok}^^i and 
{-^fe}fc=i consider Fig. ^ and define the inverse map uj = 



M) 



(z). Assume now that we already have 'I'^" 



and 



growth steps, and we want to perform the next iteration. 
To construct $(")(a-') we advance our mapping in the di- 
rection of a point z in the z-plane by adding a bump in 
the direction oi w — x^''^^^-?) the w-plane. The map 
$(")(a;) is obtained as follows: 



$(")(c^) = (0e„,A„(o.)) . (50) 

The value of 0„ is determined by 

0„ = arg[x("-^Hz)] (51) 

The magnitude of the bump A„ is determined by requir- 
ing fixed size bumps in the z-plane. This means that 



A. 



Ao 



(e' 



(52) 



We note here that it is not necessary in principle to have 
fixed size bumps in the physical domain. In fact, adap- 
tive size bumps could lead to improvements in the pre- 
cision and performance of our scheme. We consider here 
the fixed size scheme for the sake of simplicity, and we 
will show that the accuracy obtained is sufficient for our 
purposes. Iterating the scheme described above we end 
up with a conformal map that is written in terms of an 
iteration over the fundamental maps (|47|l : 



$(")(w) = 4,Q^.^^ 0...0 (/)e„,A„(w) • 



(53) 



For the sake of newcomers to the art of iterated conformal 
maps we stress that this iterative structure is abnormal, 
in the sense that the order of iterates in inverted with 
respect to standard dynamical systems. On the other 
hand the inverse mapping follows a standard iterative 
scheme 



X(")(z)=0, 



(54) 



therefore also its analytic inverse x^" (z) after n — 1 



The algorithm is then described as follows; first we 
divide the curve into segments separated by points {zi\. 
The spatial extent of each segment is taken to be approx- 
imately -s/Ao , in order to match the size of the bumps in 
the 0-plane. Without loss of generality we can take one 
of these points to be at the center of coordinates and to 
be our starting point. From the starting point we now 
advance along the shape by mapping the next point Zi on 
the curve according to the scheme described above. The 
resulting map <I>(")(z/;) is employed as 03 above. 



V. TWO- VOID MODELS: RESULTS AND 
DISCUSSION 

A. Computational aspects and Precision 

The method presented in the previous section has two 
stages which involve numerical approximation. The first 
is a numerical Fourier transform of the conformal map 
(see Eq. I34f) and the second is the truncation scheme 
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FIG. 12: The configuration of a fine with a circie. The com- 
parison of results exhibited in Fig. 1131 refers to this configu- 
ration. 
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FIG. 13: for the configuration shown in Fig. 1121 Note 
that the results are normalized by the stress intensity factor 
with no void i.e K^g = a"^ ypKa. The stress intensity factor 
was calculated near the point A (i.e. at the tip which is close 
to the void). The results were calculated for a/{d~ R) = 0.4, 
varying d/R in the range [0.1,0.7]. Isida's results (b) are 
shown in squares, and ours (a) in circles. 



(see discussion after Eq. (|35|l 'l. Aside from these two 
steps the method is analytical. In appendixIXlwe test the 
truncation scheme (with no numerical Fourier transform) 
for the case of two circles and find it to be extremely 
precise. In the next subsection we present a comparison 
of our method to another theoretical calculation for the 
straight crack and void geometry. This comparison serves 
as a testing ground for the numerical Fourier transform 
and truncation scheme combined. 



B. Stress Intensity Factor for Straight Crack and 
void 



The problem of a crack of length 2a in an infinite do- 
main, subjected to a remote uniaxial load Uyy = a°° and 
traction-free crack faces is considered as the canonical 
problem in the theory of linear elasticity fracture me- 
chanics. The tensile stress component along the tangent 
to the crack at the tip is given by [2^ 



(55) 



is known as the mode I stress intensity factor and for 
a straight crack it is given by 



(56) 



Introducing a void in front of the crack causes an increase 
in the stress intensity factor. The extent of this increase 
depends on the void's distance from the crack and its 
radius (see Fig. 112(1 . The problem of calculating the 
stress intensity factor for a straight crack and void in 
front of it has been solved using perturbational analysis 
by M. Isida 27] . Using our method one can calculate the 
stress intensity factor for the configuration in Fig. 1121 bv 
calculating the stress in a region close to the crack tip 
and fitting it to the form given in H55|) . A comparison 
of our results and those of Isida is given in Fig. We 
note that Isida's results where extracted by hand from 
a graph. Also, there are two small differences between 
our geometry and that of Isida. First, the crack is not 
strictly a branch cut but has a finite radius of curvature 
at the tip and second, in our case the inclusion deviates 
slightly from a perfect circular shape. All these factors 
together lead to an expected difference of ~ 2% between 
the two methods. We conclude that our results agree 
(to the expected precision) with those of Isida's and are 
found to be accurate even in the vicinity of the crack tip. 



C. Stress field near crack tip and void 

We now proceed to calculate the full stress field for a 
configuration of a 'rough' crack with a void. To isolate 
the effect of the void on the stress field we first calculate 
the stress for a crack without a void as before. We then 
add a void and calculate the new stress field. The hy- 
drostatic pressure and yield stress without the void are 
shown in Fig. ^1 The same quantities after the void 
insertion are shown in Fig. 1151 

In theory we could now continue the void nucleation 
process by inserting a new void in an appropriate point 
on the new yield curve (the bold curves in Fig. I15|l . It 
is obvious however that this will gain us very little. The 
new yield curve is only locally distorted by the presence 
of the void, and there is definitely no typical distance of 
2^c that could be used to create a decent two- void model. 
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FIG. 14: (color online). Panel (i); field lines of -v/jJi cf. 
Eq.Q, given in units of The bold line is the yield curve 
(i.e. \fJ2/a°° = a^/a°° = 12). The crack covers the x inter- 
val [-lO"*, 1.01 X 10*] and the bump radius is = 2. The 
loading is Mode I with ayy{oo) = o"°°. 

Panel (ii): field lines of the hydrostatic pressure, P = | Tr cr 
for the same crack as the figure above, in units of a°° . The 
bold line is the yield curve from the previous figure. A rea- 
sonable range for Pc in such a configuration is 12 < < 20, 

in accordance with ^a^<Pc<V^a-^ (see section irra . Note 
that Pc in this range ensures forward growth in the next step. 

The root of the problem is that we cannot determine 
the position of the yield curve based on a fully elastic 
solution. The correct position of the yield curve and the 
correct value of the stress on it can only come from a 
solution of the full elastic-plastic boundary value prob- 
lem. When growing the first void we assumed that the 
yield curve calculated from the solution of the fully elas- 
tic problem is close to the elasto-plastic yield curve. This 
assumption works well for the stress field of a crack 
alone, but produces physically unacceptable results for 
the stress field of a crack and void. Indeed, finite ele- 
ments calculations that take plastic fiows into account 
|16| indicate very clearly that the new yield curve is fur- 
ther removed from the first void and certainly does not 
coincide with the void boundary as it does in our calcu- 
lation. 

The unavoidable conclusion is that the doubly con- 
nected conformal calculation that is developed here is 
useful if one wants to compute the stress field of an elas- 
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FIG. 15: (color online) The same crack as in Fig. 1141 with 
the addition of a void on the yield curve. 
Panel (i) : field lines of ^/J2 in units of cr°° . The bold line is 
the yield curve with the same value of as in Fig. I14| panel 
(i) . The dotted line is the yield curve from Fig. 1141 panel 
(i) i.e. the yield curve for the same crack with no void. The 
addition of the void creates a very local perturbation of the 
yield curve. 

Panel (ii): field lines of the hydrostatic pressure, P = |Trcr 
in units of a°° . The bold line is the yield curve from the figure 
above. Again the perturbation of the field lines with respect 
to Fig. 1141 panel (ii) is very localized. 

tic material in which a hole was inserted in the vicinity 
of a crack. It cannot be used however to develop an 
approximate method of taking into account the plastic 
yields that result in a successive appearance of two voids. 
The first void can be inserted on the basis of elastic cal- 
culations, but the second void cannot be added without 
a considerably improved consideration of the plastic dy- 
namics. As long as the analytic aspects of plastic dy- 
namics are not elucidated better, the one- void model is 
proposed as the best available approach to roughening 
via growth with plastic deformations. 



VI. SUMMARY AND CONCLUSIONS 

The crack propagation model presented above, repro- 
duces the experimental appearance of self affine crack 
rupture lines with an anomalous Hurst exponent. The 
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long range correlations are created by the stress field 
which satisfies boundary conditions on the crack's in- 
terfaces. The ability to solve the elastic boundary value 
problem for an arbitrary crack is the basic building block 
of our theory which enables us to capture these correla- 
tions. 

To gain a sizeable scaling range with anomalous ex- 
ponent our model needs to employ a sufficiently wide 
probability distribution function in the angle 6 around 
the tip. While we could not discern a strong dependence 
of the numerical value of the scaling exponent on the pdf, 
we do observe a significant dependence of the extent of 
the scaling region. Whenever we found a scaling range 
the numerical value of the exponent was in the range of 
C = 0.66 ±0.03. 

We have also presented a general method for the cal- 
culation of the elastic stress field surrounding a doubly 
connected region, like a crack and void at its tip, us- 
ing conformal maps. This method, although not useful 
for the creation of a two-void model of crack growth, is 
quite general and can be used for solving different phys- 
ical problems in doubly connected domains. 
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APPENDIX A: TWO CIRCULAR HOLES 

In this section the specific case of two circular holes is 
solved using the formalism of section IIVI For the case 
of two equally sized circles there exists an analytical so- 



lution [2J|. In this solution bipolar coordinates are used 
and a series expansion for the stress components is given 
in which all the coefficients have a given closed form. We 
take advantage of the bipolar method to check the pre- 
cision of our doubly connected method in the two circle 
limit. 

The conformal map for two circular holes is. 



1 



LO — a 

a e 3fi and < a < 1. 



(Al) 



maps the annulus, i.e. p < < 1, onto the ex- 
terior of the unit circle and an additional circle on the 
right, see Fig. ^] Note that <I>(a;) is the first auxiliary 
map described in subsection IIV CI Since for this case the 
conformal map has such a simple form, one need not use 
a numerical FT to obtain its Fourier coefficients, there- 
fore the only approximation left in our method is the 




*(tB) 
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FIG. 16: An example of a conformal mapping of tlie annulus 
onto two circular holes. In this example x\ — 1.1 and X2 = 
3.1. Note that the unit circle is mapped onto itself, while the 
p-circle in the tj-plane is mapped onto the circle on the right 
in the z-plane. 



truncation of the functions and V'(^) ^^s explained 

in subsection IIV Bl As a result the two circle case allows 
us to isolate and estimate the error associated with the 
truncation approximation by comparing with the exact 
results of the bipolar method. 

Specializing Eqs. ^^-^^ for $(w) given in Eq (|X2)) 
we obtain, 



oo -. / -. N 3 oo oo 

+ E E n^ne—'+ Y: i'nS-^D,, (A2) 

n— — OO n— — oo n— — oo 

and 

^f^^ + ^)+ E E n^^p-h—^+ ± ^„p«e-" = i^.. (A3) 

2 \ pe — a p — as I ^-^ \ — pe — a ^-^ ^-^ 

^ ^ n— — oo n— — oo n— — oo 

I 

As mentioned above, by expanding these equations in this method we calculated the stress field surrounding 
powers of e one can get a well-posed set of equations two equally sized circles see for example Fig. ^] We 
which can be solved by the truncation procedure. Using then solved the same problem using the fully analytical 



cr°° / ae — 1 e — a \ 
2 \ e — a ae — 1 J 
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bipolar method and found that both methods give the 
same resuhs to a very high degree of accuracy, i.e. in 
all space the difference, Act, between the two methods 
satisfies 

^ - 10-- . (A4) 



This holds in particular for area between the voids which 
might be expected to give convergence problems. We 
conclude that the truncation method proves itself very 
accurate in this limiting case. 
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